////////////////////////////TCF/////////////////////////////
// TPC Coherent Fit (TCF)                                 //
// Author: Xianguo Lu, lu.xianguo@gmail.com               //
// Reference: Xianguo Lu, "Exploring the performance      //
// limits of the ALICE Time Projection Chamber and        //
// Transition Radiation Detector for measuring identified //
// hadron production at the LHC", dissertation at the     //
// University of Heidelberg (2013).                       //
////////////////////////////////////////////////////////////

#ifndef _MN2D_H_
#define _MN2D_H_

#include "dEdxGausN.h"

class MN2D
{
  //==============================================================================
  //                                    controls
  //==============================================================================
 public:

 private:
  static Double_t MaxIter(){return 1e8;}

  static Int_t fkParaResolution;
  static Int_t fNPolint; //must be an even number
  static Double_t fTolerance; 
  static Int_t fRegOpt; 
  static Bool_t fkFixMean;
  static Double_t fElectronHighPTrigger;
  static Double_t fElectronHighPFrac;
  static Double_t fRegWeight;
  static Double_t fCpIni;
  static Bool_t fkUseTrueYield;

  //==============================================================================
  //                                    initialization and destruction
  //==============================================================================
 public:
  static void Ini(TH2D *hpreID, TH2D * horiginal,  const Int_t nb, const Int_t binx[], const Double_t bbpar[], const Double_t binpar[], const Double_t realRes);
  static void Clear();
  
 private:
   
  //==============================================================================
  //                                    functions
  //==============================================================================
 public:

 private:
  static Double_t IDdEdxM1Q0(const Double_t *par, const Double_t mean);
  static Double_t IDdEdxM2Q0(const Double_t *par, const Double_t mean);
  static Double_t IDdEdxM3Q0(const Double_t *par, const Double_t mean);
 
  //==============================================================================
  //                                    parameters
  //==============================================================================
 public:
  static Int_t NparBin(){return dEdxGausN::Npar();}

 private:
  static Int_t Npar(){return NparBuffer()+fNbin*NparBin();}
  static Int_t Nfree(){return Npar()-fNfix;}
  static Int_t NparBuffer();
  
  static void IniResPar(Double_t * respar, const Double_t realRes);
  static Int_t NparRes();
  static Double_t GetResolutionParametrization(const Double_t mean, const Double_t * par);
  
  static void FixOverflow();
  static void SetIfixed();
  static Int_t SetIfixAllRes(Bool_t pfix[], const Bool_t kfix);
  static void SetLowHigh();
  static void FixHighPElectron();

  static const Double_t * GetBinPar(const Int_t ibin, const Double_t xx, Double_t *localpar, Int_t *outorder, const Bool_t kprint=0);
  static void SetPhysics(const TMatrixDSym * fullcov, const Bool_t khesse);

  static Int_t fNfix;
  static Bool_t *fIfixed;
  static Double_t *fLows;
  static Double_t *fHighs;
  static Double_t *fPars;
  static Double_t *fLowErrs;
  static Double_t *fHighErrs;
  static Double_t *fAs;
  static Double_t *fErrAs;
  static Double_t *fRatio;
  static Double_t *fErrRatio;
  static Double_t *fACov;
  //==============================================================================
  //                                    output
  //==============================================================================
 public:
  static void Draw2d(const TString tag);

 private:
  static Char_t* Status(Char_t * info, const Bool_t kend=kTRUE);
  static void Print(const TString tag, const Bool_t kprintBinPar=kTRUE, const Bool_t kprintBinErr=kTRUE, const Bool_t ktofile=kTRUE);
  static void PrintPhysics(FILE *pout, const TString ss, const Double_t * var);
  static void Save(const TString tag);
  static void SavePhysics(TTreeSRedirector * stream, const TString tn, const Int_t ibin);
  static FILE * fLog;

  //==============================================================================
  //                                    fit
  //==============================================================================
 public:
  //
  static void Fit(const Int_t startnbeta=0);

 private:
  static Double_t ScanShift(const Double_t standardpars[], Int_t &beta, const Int_t itarget, const Double_t p1, const Double_t p2, const Double_t dp);
  static void ScaleAll(const Double_t standardpars[], Int_t &beta, const Double_t scale);
  static void DoShift(const Double_t standardpars[], Int_t &beta, const Int_t itarget, const Double_t corr);
  static void FinalFit(Int_t &nbeta);

  static Double_t LinearFit(const Int_t np, const Double_t xarr[], const Double_t yarr[], const Double_t eys[], const Double_t xtarget, const Int_t order);
  //
  static Double_t RegularizationWeight(const Double_t ncount);
  static Double_t RegMaxdEdx(){return 110;}
  //
  static void UpdateFractionWithPreID(Double_t &p0y, const Double_t tpcall, const Double_t prepar, const Double_t preall);
  //
  static Double_t GetFCNChisquare(Double_t pars[]);
  //
  static void MinFCNChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *pars, Int_t iflag);
  //  
  static Int_t DoMigrad(const Int_t nbeta, const Bool_t klong);
  //
  static void BetaFit(Int_t & nbeta, const Bool_t kHESSE);
  //
  static Double_t SetCrossing();

  static Int_t NDF(){return fNdata - Nfree();}
  static Double_t ReducedChi2(){ if(!NDF()) exit(1); return fStatChi2/NDF();}
  static Double_t myFCN(){ return fStatChi2 + fBeta*fPenalty; }

  static TH2D *fHoriginal;
  static TH2D *fHFitData;
  static TH2D *fHprefilter;
  static TH2D *fHrest;
  static TH1D *fHFitYield;
  static TH1D *fHTrueYield;
  static TH1D *fHrestyield;

  static Int_t fNbin;
  static Int_t fNYbins;
  static Int_t *fBinX;
  static Double_t *fXCenter;
  static Double_t *fFitYieldArray;
  static Double_t *fTrueYieldArray;

  static Double_t fMinXing;
  static Double_t fMaxXing;

  static Bool_t fkPreID;
  static Double_t *fPrePar;
  static Double_t *fPreAll;

  static Double_t fYbw;
  static Double_t *fYCenter;
  static Int_t fNdata;

  static Double_t fStatChi2;
  static Double_t fBeta;
  static Double_t fPenalty;

  static TStopwatch * fStopwatch;
  static TMinuit *fMinuit;
  static Int_t fkErr;

  //==============================================================================
  //                                    end
  //==============================================================================
};

#endif
